Peer Community Journal

In search of recurrent patterns of postzygotic hybrid incompatibility, we investigated the Holarctic 'Old World Swallowtail' (Papilio machaon Linnaeus) butterfly complex, whose members can easily be crossed in the laboratory. How many species this model system comprises remains unclear, as taxa with highly distinctive larvae but uncertain status come into contact with authentic P. machaon subspecies in Southern parts of the Palearctic region. By determining mitochondrial and ITS2 haplotypes within and away from contact zones, we found that these neighboring populations do exchange genes, as expected from F1 hybrids being generally fertile in our laboratory crosses. Nevertheless, recurrent instances of dysregulation of diapause were uncovered in hybrid progenies. In keeping with Haldane's Rule, pupae of the heterogametic (female) sex were either unable to enter diapause or, in reciprocal crosses, unable to resume development after having initiated diapause, whereas F1 males experienced normal, photoperiod-regulated diapause, but passed on abnormal diapause regulation to part of their female offspring when backcrossed. Comparing male and female pupal weights in hybrid progenies provides estimates of parental incompatibility that allow to rank taxa and predict quantitatively the outcome of additional crosses, as expected if the same regulatory system were involved. A survey of the entomological literature confirms that diapausing pupae that cannot resume development ('perpetual nymphs') are a recurrent feature of interspecific lepidopteran crosses. Moreover, of the two parent species of a perpetual nymph, the paternal one generally has fewer broods per year. These observations are discussed in the light of models of the speciation process.

We describe a chromosome-level, haplotype-resolved genome assembly from a male European moose (Alces alces alces). The assembly comprises two pseudo-haplotypes of 3,148 Mb and 3,112 Mb, with sex chromosomes in haplotype one, and 33 autosomes in each haplotype (68 in total). Assembly completeness is high (BUSCO 98.3% and 95.7%), with 21,496 and 20,498 annotated protein-coding genes for haplotypes one and two, respectively. This genome assembly is the most complete so far generated for European moose.

In trypanosmatids genes, transcribed by RNA polymerase II do not have canonical promoters and are organized into directional gene clusters that mature into monocistronic transcripts by a co-transcriptional process known as trans-splicing. Even though gene expression is regulated mainly post-transcriptionally, it is currently understood that chromatin and epigenetics are also involved in this regulation. In eukaryotes, specific signals are normally required for the occurrence of an appropriate transcription initiation. Among them, trimethylation of histone H3 in lysine 4 is the most conserved signal normally detected at transcription start sites of actively transcribed genes. Unlike many model organisms, trypanosomes do not have defined promoters. Instead, transcription initiates in a bidirectional manner from dispersed regions coincident with divergent strand switch regions located between directional gene clusters (DGCs). In T. cruzi, H3K4me3 was observed at the origins of transcription coincident with divergent strand switch regions (dSSRs) in epimastigotes, but it has not been mapped throughout the whole genome at base-pair resolution or in other life stages so far. Here, we set up the CUT&RUN technique for T. cruzi epimastigotes and trypomastigotes. Consistent with a predominant post-transcriptional regulation along the life cycle, we did not find significant differences between life stages. We corroborated that H3K4me3 is enriched at dSSR adjacent to actively expressed DGCs. Moreover, we noticed that this histone mark exhibits different patterns that correlate with the genomic span of the transcription initiation regions and with transcriptional activity. Furthermore, we unveiled that the most actively transcribed DGCs are associated with shorter dSSRs and are located within the core compartment of the genome displaying a more accessible chromatin.

BackgroundLyme borreliosis is a common tick-borne disease in Europe caused by spirochetes of the Borrelia burgdorferi sensu lato complex, including Borrelia afzelii, which is maintained in nature through interactions with rodent reservoir hosts. These spirochetes have evolved several surface proteins to manipulate rodent host immunity, some of which remain polymorphic in Borrelia populations. Among these proteins, OspE, which binds the host complement-regulating factor CFH to evade destruction by complement, is one of the most variable. Yet, what evolutionary forces maintain this polymorphism is not well understood. Motivated by a recent discovery of CFH polymorphism in the bank vole (Clethrionomys glareolus), the main reservoir host of B. afzelii, we hypothesized that the polymorphism is maintained by host-parasite coevolution involving specific associations between host and parasite genetic variants. MethodsWe analyzed associations between bank vole CFH alleles and B. afzelii OspE variants across three datasets sampled in Poland. Selection acting on OspE was evaluated using omegaMap. Host-pathogen genotype associations were tested using partial redundancy analysis (RDA), and co-structure was assessed using co-correspondence analysis (CoCA). ResultsWe found that OspE evolves under positive selection, however, we found no evidence for an association between OspE and host CFH variants at the individual level based on RDA or at the population level based on CoCA. ConclusionsDespite evidence of positive selection acting on OspE, we found no support for specific genetic matching between B. afzelii and its bank vole host at the CFH-OspE interface. These results suggest that the evolution of CFH and OspE may be shaped by broader selective pressures, potentially including interactions with multiple host species.

The migratory endoparasitic pinewood nematode (PWN), Bursaphelenchus xylophilus, is the causal agent of pine wilt disease, causing significant economic and ecological losses in conifer forest ecosystems in Europe and Asia. Understanding the molecular mechanisms regulating PWN parasitism-related genes may lead to new sustainable solutions for control. Based on previous PWN transcriptomic datasets from the pre-parasitic and parasitic stages and from the pharyngeal gland cells (GC), an in silico analysis was performed to identify transcription factors (TF) highly expressed in the GC. Seven candidates TF genes were selected, and their spatial expression validated by in situ hybridisation. From those, two GC-expressed TFs, BXY_079 and BXY_022, each encoding zinc finger domains, were successfully knocked down by RNA interference. Transcriptomic data from silenced BXY_079 and BXY_022 TFs, analysed with existing life cycle specific transcriptomic data, showed that both TFs control genes expressed at similar times, by repressing male-related genes while activating genes expressed during the J3 and D3 stages, yet each represents the extreme of the others minor function. In addition to these common roles, BXY_079 also activates parasitism-related genes in the J2 stage. These BXY_079-activated parasitism-related genes predominantly encode proteins with lytic functions, including secreted peptidases and glycoside hydrolases. Consistent with their proposed role in parasitism, these genes are highly expressed during the parasitic juvenile stages and are likely involved in nematode feeding, tissue penetration, and migration within the host. In contrast, BXY_022 also represses the expression of several genes related to the reproduction system, such as major sperm proteins and cytosolic motility proteins, particularly in the adult male stage. Taken together, both dual-functional TFs work together, non-redundantly, to regulate gene expression across the life cycle, while each is additionally specialised to regulate diverse and distinct gene sets: ranging from genes implicated in lytic parasitic functions to sexual dimorphism.

O_LIStudying animal movement at the population scale requires a stable, modern software substrate. Within R, the legacy move package supplied that substrate for over a decade, but its sp/rgeos backbone has been retired. The successor package move2 deliberately confined its scope to the data class and core movebank API functions. C_LIO_LIThe analytical machinery of move, namely dynamic Brownian-bridge utilisation distributions, the directional bivariate-Gaussian variant, corridor segmentation, and along-track thinning, was left to port to the modern sf/terra stack. C_LIO_LIWe present move2utils, an R package that completes and complements that transition. move2utils provides move2-native ports of the move analytical functions, preserves the original C kernels where they exist, and replaces the deprecated spatial scaffolding around them. It additionally ports some of the legacy R-based code to faster C kernels to improve computational speed. move2utils also exposes novel outlier-detection methodology described in detail in a companion paper. C_LIO_LIThe package is open-source (GPL [≤] 3), is developed on the MPCDF GitLab and mirrored on GitHub for public installation, and ships with vignettes and a CI-tested check suite. We illustrate it with a worked example on real tracking data and synthetic datasets. C_LI

Arctic and Boreal raptor communities will continue to be affected by borealization and other climate change related processes, providing a challenge for ecologists predicting future sates. However, by using community assembly theory and species traits, future communities may be predictable. In this study, we analyzed variation in reproduction traits as a consequence of diet specialization for 29 raptors, 2 skuas and 3 corvids. We assessed and implemented foraging traits for specialists and generalists into predator-prey models from which successful invasion conditions were derived. Specialist raptors produced larger clutch sizes, had a higher proportion of fledged per clutch and also expressed more variation compared to generalist raptors. These results suggest a relationship between diet specialization and reproductive traits which was also observed within phylogenetic orders. Specialist owls (Strigiformes) produced higher clutch sizes with a larger clutch range compared to generalist owls. The same pattern was observed for falcons (Falconiformes). No clear difference in reproduction was observed for specialist and generalist hawks, kites and eagles (Accipitriformes). Corvids expressed clutch sizes similar to that of specialist raptors while having the lowest proportion of fledged per clutch. Differences in foraging traits between specialists and generalists could be distinguished using functional response curves. A predator-prey model parameterized with foraging trait data showed that a generalist can coexist with a resident specialist if it has access to prey unavailable to the resident specialist. Otherwise, the native specialist outcompetes the invading generalist due to foraging efficiency. The combined empirical and theoretical findings in this study show how diet specialization affects both reproduction and the potential invasion success of raptors.

While emerging diseases pose a global threat to amphibians, the dynamics of understudied vector-borne blood pathogens remain poorly understood. Pathogen occurrence is driven by a combination of environmental, ecological, and phylogenetic factors, yet how these drivers shape blood pathogen communities in tropical amphibians is largely unknown. In this study, we used molecular screening and phylogenetic linear models (PGLMMs) to evaluate how climate and ecomorphology influence the incidence of three blood pathogen groups (Trypanosomatidae, Hepatozoon, and Rickettsia) in wild anurans from a protected area in the Brazilian Atlantic Forest. Among 93 individuals sampled, over 93% were infected with at least one pathogen. Trypanosomatidae was the most common (76.3%), followed by Rickettsia (69.9%) and Hepatozoon (16.1%). Pathogen responses to temperature were contrasting: Hepatozoon occurrence increased in warmer periods, while Trypanosomatidae declined. Furthermore, rheophilic species showed a lower probability of Rickettsia infection, providing the first evidence that microhabitat use influences blood pathogen dynamics in amphibians. Our findings demonstrate that hemoparasites prevalence is driven by a multifaceted interplay of variables, highlighting that conservation strategies must account for these pathogen-specific responses to habitat use and environmental change, even within protected areas.

The climatic oscillations of the Quaternary have likely affected the demographic history of many species, and PSMC (Pairwise Sequentially Markovian Coalescent) has been widely used to investigate these histories. However, it is increasingly acknowledged that PSMC trajectories are difficult to interpret. First, they are influenced by connectivity changes, even without population size changes. Second, most PSMC curves exhibit a few humps when tens of cycles occurred during the Pleistocene. Finally, responses to ancient habitat change have been shown to be species-specific. To address these issues, we simulated structured populations where connectivity (or population size and connectivity) varied according to successive interglacial and glacial periods during the last 2.6 million years. We computed the IICR (Inverse Instantaneous Coalescence Rate), the function that PSMC estimates, and ran PSMC. We further varied the generation length and assumed that some species were positively or negatively affected by glacials. We found that the IICR carries information regarding the demographic oscillations, but that PSMC fails to recover it for times older than 300 ky. For the last 200 ky, PSMC was often able to reproduce qualitatively the demographic oscillations. We also tested SNIF (Structured Non-stationary Inferential Framework), which produced good results using the IICR curve as an input but not when using the PSMC curve. Altogether, our study suggests that the humps older than 300 ky in PSMC histories are unlikely to represent trends of population size or connectivity. However, improving the estimation of the IICR could potentially help reconstruct some of these past demographic changes.

BackgroundThe invasive mosquito Aedes koreicus has established populations in several European countries during the past decade, raising increasing public health concerns due to its potential role as a vector of pathogens. While species identification is primarily based on morphological characters, Ae. koreicus exhibits distinct morphological variants originating from mainland Korea and Jeju Island, which complicates surveillance and may lead to misidentification, particularly in regions where closely related species co-occur. To date, the genetic basis and population-level relevance of these morphological forms in Europe remain poorly understood. MethodsWe investigated the co-occurrence of two morphological forms of Ae. koreicus in Hungary, representing the first confirmed European location where both forms were detected sympatrically and even within the same breeding sites. Adult mosquitoes were morphologically characterised using diagnostic traits, and individuals representing both morphotypes were subjected to comprehensive genetic analyses. We sequenced multiple mitochondrial markers (COX1, COX2, COX3, ATP6, ND1, ND3) and the nuclear ITS2 region using Oxford Nanopore long-read sequencing. Phylogenetic reconstructions and haplotype network analyses were applied to assess genetic differentiation between morphotypes and to compare them with Aedes japonicus as a closely related reference species. ResultsAcross all analysed mitochondrial and nuclear markers, no genetic differentiation was detected between specimens identified as the "mainland" or "Jeju-do" morphological forms of Ae. koreicus. Phylogenetic and haplotype network analyses consistently grouped individuals independently of their morphotype, indicating a shared genetic background at the population level. In contrast, Ae. japonicus formed a clearly distinct genetic lineage, confirming the robustness of the applied markers for interspecific discrimination. ConclusionsOur results demonstrate that the observed morphological variability in Ae. koreicus populations are not underpinned by detectable genetic differentiation using commonly applied mitochondrial and genomic markers. These findings highlight the limitations of using morphology alone to infer population origin or structure and emphasise the need for heightened awareness of intraspecific variability in routine surveillance. Accurate morphological identification remains critical, particularly in citizen science-based monitoring programmes and AI-based automatization platforms, such as Mosquito Alert, to avoid confusion with morphologically similar invasive species. Further studies integrating genomic and ecological approaches are required to elucidate the mechanisms underlying morphological variability in this emerging vector species.

Phylogenetic reconstruction is a fundamental problem in comparative genomics. As a theoretical problem in rearrangement studies, this has been modelled as the Small Parsimony Problem (SPP), in which ancestral genome structures have to be determined minimizing the number of rearrangement events occurring throughout the phylogeny. This problem is of significant interest in microbial and cancer genomics, due to the prevalence and clinical importance of rearrangement events. Genome structures in this problem are expressed as sequences of markers, which are themselves oriented sequence features (such as genes) that abstract from non-structural variations. Recent research has focused on the problem under the natural genomes model, in which arbitrary variations in copy number of markers are allowed. Natural genomes are often studied under the DCJ-indel model, a model which has already been successfully applied to plasmid data. There also exist ILP solutions to a variant of the Small Parsimony Problem under the DCJ-indel model. However, these solutions are limited in their applicability, as they make some critical simplifications for tractability purposes: ancestral marker frequencies and precomputed putative ancestral adjancencies, with their predicted likelihoods, are assumed as input. This creates multiple problems from both a theoretical and practical perspective. Firstly, this simplification means that not the full state space is searched for a solution, but rather only the subset of genomes with the precomputed putative adjacencies, meaning an optimal solution to the exact SPP is not guaranteed. Secondly, marker frequencies are given externally, without any theoretical guarantees. Thirdly, the method used to precompute adjacencies relies on gene trees, which requires the use of genes as markers, when gene annotation is often unreliable, especially in regions with a lot of rearrangement. Additionally, this restricts the applicability of the approach to sets of genomes that are both divergent and large enough to be able to produce informative gene trees. This is, for example, rarely the case for plasmids, where nucleotide mutations are rarer than rearrangements and genomes are small. Hence, we revisit the problem to solve the exact SPP by introducing a cost to indel operations, which allows us to compute ranges of marker frequencies and derive theoretical results, that allow us to reduce the solution space that the ILP searches without sacrificing optimality. We show that this makes the problem tractable for the case of small and recently related genomes, first on simulated genomes, and then on a set of pathogenic plasmids which represent a realistic use case for the method.

Predators are often considered regulators of disease in prey populations, a concept central to the "healthy herd hypothesis". This hypothesis suggests that by preferentially removing infected individuals, predators can reduce parasite prevalence. However, predators may also act as disease vectors, facilitating the spread of parasites. We investigated whether stickleback fish (Gasterosteus aculeatus) can act as vectors for the transmission of the obligate bacterial parasite Pasteuria ramosa to its Daphnia host, a widespread freshwater zooplanktor. We fed infected D. magna to sticklebacks, and subsequently analysed faecal samples for the presence, viability, and infectivity of parasite transmission stages (= spores). We recovered approximately 60% of the consumed spores from fish faeces and these spores did not suffer from reduced infectivity to D. magna. Additionally, spores associated with sloppy feeding did not reduce infection rates. Thus, consumption of infected hosts by fish does not eliminate the parasite, but in contrary, may contribute to the spread and persistence of P. ramosa in natural populations, potentially influencing parasite dynamics in natural freshwater ecosystems.

Lumpy skin disease (LSD) is an emerging cattle disease caused by lumpy skin disease virus (LSDV), with major impacts on the industry, being classified as a Category A disease. Although it was historically confined to Africa, LSD has expanded into the Middle East, Asia and Europe. Here, we report two LSDV genomes from the first outbreak detected in Catalonia, Spain, in October 2025. The genomes were assembled from high-throughput sequencing data generated from two homogenized skin nodules. Comparative phylogenetic analyses were performed using all available complete LSDV genomes and rpo30 gene sequences. These analyses placed the LSDV isolates detected in Catalonia within clade 1.2, closely related to the isolates recently reported in Sardinia, Italy. Our findings also support a connection between recent south-western Europe and central African strains, possibly through northern Africa, and highlight the need for more complete genomes to clarify the origin and connections among recent LSDV outbreaks.

Body size, through its links to various physiological traits, has often been hypothesized to influence evolutionary rates. Negative body size-rate correlations have been reported in the morphological or molecular evolution of several extant vertebrate groups, including mammals, birds, reptiles, and teleost fishes. In this study, we estimated body masses for 89 species of plesiosaurs, a clade of Mesozoic aquatic reptiles, and found that their body size evolution conforms to a three-regime Ornstein-Uhlenbeck process, indicative of constrained evolution. Rates of morphological evolution, inferred using the skyline fossilized birth-death process and the variable-rates model, show minimal support for a correlation with body size in this clade. Our results thus serve as a counterexample, suggesting that the negative body size-rate relationship is not a universal vertebrate pattern, but rather a trend restricted to certain lineages.

Age-dependent survival is central to understanding population dynamics and life-history evolution. We analysed carcass weight and age-at-harvest data from 6022 red foxes (Vulpes vulpes) collected across Sweden between 1967 and 1971 to evaluate latitudinal effects on body mass and age-dependent survival. Carcass weights decreased from south to north in both adults and sub-adults, contrary to Bergmann's rule, with southern foxes weighing approximately 1.27 times more than northern foxes. The latitudinal weight gradient exceeded the sex difference in both age classes, and no sex x region interaction was detected. The decrease in weight with latitude is consistent with reduced prey availability and harsher winter conditions in the north, which limit growth and body size during development. Using a Bayesian age-at-harvest model with region-specific population growth rates (lambda), we estimated age-dependent survival probabilities for four latitudinal regions and both sexes. Despite the strong latitudinal gradient in weight, survival did not show a corresponding pattern - regional differences were uncertain, with all credible intervals spanning zero. Regional population growth rates were consistent with slight decline in the north and near-stability in the south-central region, which suggests that body condition and population dynamics are coupled at the regional scale despite no survival gradient. The decoupling of body condition and survival across regions suggests that mortality patterns are similar across the latitudinal gradient. We discuss these patterns in terms of latitudinal productivity gradients, prey availability, and life-history trade-offs in a widely distributed carnivore.

Free-roaming dog populations pose challenges to public health and animal welfare, and their management and protection is legally mandated in many countries to control zoonotic diseases such as rabies, reduce bites, mitigate fear of aggressive dogs, and preserve animal welfare. Although international organisations recommend holistic dog population management with surgical sterilisation as a key component, stakeholders report a lack of evidence-based guidance on effective and cost-efficient strategies. To address this, we developed a deterministic model of a closed free-roaming dog (FRD) population, including dogs dependent on and independent of humans, parameterised using data from southern India. We evaluated a wide range of sterilisation strategies varying duration, interval, frequency, intensity, targeting different population subsets, and fixed or responsive approaches. Strategies were assessed based on their effects on independent FRD population size and implementation costs. Female sterilisation coverage (FSC) strongly determined population reduction and cost-effectiveness. Repeated, short sterilisation periods with appropriately timed intervals and high early intensity best balanced population reduction and cost-efficiency. Excluding dependent dogs (cared for by humans) severely constrained population reduction, while male sterilisation enhanced impact at little extra cost, assuming direct effects on pregnancy rates. Costs and population impacts can be visualised on an interactive web application (https://field.shinyapps.io/DogPopSimApp), allowing practitioners and policymakers to explore sterilisation strategy outcomes. Overall, cost-effective population reduction required sustained commitment and coverage-driven strategies that balance welfare, population reduction and the operational needs of dog-mediated rabies control.

Various hypotheses have been proposed to explain why some species persist or even flourish in urban areas. Yet, despite its central role in determining when and where animals encounter resources, disturbance, and risk, movement behaviour remains an overlooked mechanism of urban success. In urban areas, human activities are strongly periodic, i.e. predictable in space and time. This may favour species able to adjust their behaviour to predictable cycles of resources and risks in space and time. Here, we tested this hypothesis and tracked movement behaviour along an urbanisation gradient in three mammal species with different urban success: red fox (Vulpes vulpes), an urban dweller; raccoon (Procyon lotor), an invasive urban dweller; and wild boar (Sus scrofa), an urban utiliser. We analysed periodicity in movement behaviour and investigated whether increasing urbanisation is associated with periodic reorganisation of activity timing, space use, and further analysed alterations in habitat selection along the urbanisation gradient. Our results show that foxes aligned their movement behaviour with human activity, having stronger day-night contrasts and more repeatable space use than their rural counterparts. Urban raccoons showed a contrasting strategy; they were more active during the day, without changes in their movement routines under increasing urbanisation, suggesting a flexible strategy that explains their urban success. In contrast, wild boars reduced routine movement behaviours with increasing urbanisation, consistent with their occurrence in less predictable suburban environments and avoidance of city centres. In summary, our results suggest that movement behaviour may be a key mechanism enabling animals to persist in cities, revealing distinct behavioural strategies for coping with urban environments.

The study of historical documents from the original excavation led by J. M. Cruxent, of museum collections, as well as results of fieldwork, are combined to provide a revised assessment of the archeologically site of Cucuruchu, which together with those of Taima-Taima and Muaco, document megafaunal lithic remains associated from the Late Pleistocene in northern Venezuela. Six vertebrate taxa are reported for Cucuruchu: five megaherbivores (a giant sloth, Eremotherium laurillardi; a glyptodontid, Glyptotherium cf. G. cylindricum; a macrauchenid, cf. Xenorhinotherium bahiense; a proboscidean, Notiomastodon platensis; and an equid, Equus sp.), and a turtle (Kinosternon sp.). Two fragments of El Jobo projectiles were found in situ lying among megaherbivores bones, while another unidentified projectile fragment was also found in the same layer. Although no direct evidence of human exploitation on these megaherbivores remains has been found, the presence of El Jobo projectiles in the undisturbed fossil-bearing layer some of which are associated with fossil remains, could be a suggestive of such interaction. Two new radiocarbon dates indicate ages of at least ca. 15.3-16.0 kyr cal BP for the fossil-bearing layer, contrasting with the mid to late Holocene age originally reported. Cucuruchu is a window to the past and supports the presence of hunter-gatherers of the El Jobo culture in the region at the end of the Pleistocene.

High concentrations of steroidal hormone compounds are a growing source of concern for environmental pollution in aquatic ecosystems. In this study, we examine the effects of two estrogenic compounds (estriol and 17-ethinylestradiol) on fitness traits in the aquatic microcrustacean, Daphnia magna, a key bioindicator species for toxicology studies. The impacts were compared of two forms representing a natural and synthetic estrogenic compound. Growth and reproduction traits were assayed by exposing Daphnia to each estrogen type at four concentrations reflecting potential environmental exposure conditions up to acute toxicity levels (ranging from 0.1 - 50 {micro}g/L). Assaying the effects at a variety of concentrations is important given that it is known that hormone exposures can often result in non-monotonic responses. Both forms of estrogen impact a subset of the traits assessed, in some cases leading to beneficial changes and others causing harm. Estriol, the naturally-occurring estrogen, and EE2, the synthetic version, at high doses shift fitness traits in opposite directions such as adult growth rate as do at low doses for fecundity. In conclusion, our results support the need to assay a wide array of traits using multiple forms of steroidal hormones at a range of doses in order to assess non-monotonic patterns and their impact on an organismal fitness. In particular, assays that extend beyond the conventional measurements of lethality during acute exposure windows will be essential for understanding the impact of increased levels of hormone pollution on aquatic organisms and ecosystem health.

Comparative transcriptomic datasets are increasingly used to investigate the molecular basis of phenotypic diversification across species. However, finding genes that are differentially expressed (DE) between lineages remains challenging, for two main reasons. First, the random evolutionary drift can blur the signal left by lineage-specific shifts in mean expression, and induces phylogenetic correlations that, if ignored, can widely inflate the False Discovery Rate (FDR), i.e., the amount of spuriously detected genes. Second, DE analysis from RNA-Seq data involves multiple testing on many genes for a small number of individual measurements with high noise, and requires dedicated statistical tools. Traditional DE tools, such as limma, and classical Phylogenetic Comparative Methods (PCMs), such as the Expression Variance and Evolution (EVE) model, are both designed to tackle one of these two challenges alone, but both fail in the context of inter-species RNA-Seq data. In this work, we present phyloDE, a new tool for inter-species DE, that aims at taking the best from both approaches. On simulations based on a recently published four-species rodent dataset, we show that, contrary to other methods, phyloDE correctly controls the FDR in all settings, while keeping a reasonable power. When reanalyzing the empirical dataset, phyloDE discovers more DE genes that exhibit consistent changes in their cis-regulatory landscape compared to EVE in all the experimental settings. The method is implemented in R, with an interface inheriting from limma.